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O ' ABSTRACT 

\ We introduce STECKMAP (STEllar Content and Kinematics via Maximum A Pos- 

teriori), a method to recover the kinematical properties of a galaxy simultaneously 
" ' \ with its stellar content from integrated light spectra. It is an extension of STECMAP 

, (Ocvirk et al. 2005) to the general case where the velocity distribution of the under- 

' lying stars is also unknown. The reconstructions of the stellar age distribution, the 

\ age-metallicity relation, and the Line-Of-Sight Velocity Distribution (LOSVD) are 

■ all non-parametric, i.e. no specific shape is assumed. The only a propri we use are 

' positivity and the requirement that the solution is smooth enough. The smoothness 

\f~^ , parameter can be set by generalized cross validation according to the level of noise in 

C3 ' the data in order to avoid overinterpretation. 

We use single stellar populations (SSP) from Pegase-HR {R = 10 000, AA = 
pi • 4 000 — 6 800 A, Le Borgne et al. 2004) to test the method through realistic simulations. 

I ' Non-Gaussianities in LOSVDs are reliably recovered with SNR as low as 20 per 0.2 

O ^ A pixel. It turns out that the recovery of the stellar content is not degraded by the 

■ ■ simultaneous recovery of the kinematic distribution, so that the resolution in age 

and error estimates given in Ocvirk et al. (2005) remain appropriate when used with 
STECKMAP. 

We also explore the case of age-dependent kinematics (i.e. when each stellar com- 
ponent has its own LOSVD). We separate the bulge and disk components of an 
idealized simplified spiral galaxy in integrated light from high quality pseudo data 
' (SNR=100 per pixel, R=10 000), and constrain the kinematics (mean projected veloc- 

ity, projected velocity dispersion) and age of both components. 

Key words: methods: data analysis, statistical, non parametric inversion, galaxies: 
kinematics, stellar content, formation, evolution 



C/3 



1 INTRODUCTION 

For decades now, the spectral indices from the Lick group 
have been used to study the properties of stellar popula- 
tions (Faber et al. 1985; Worthey 1994; Trager et al. 1998). 
Since the profile and depth of the lines involved in these 
spectral indices are afltected by the Line Of Sight Velocity 
Distribution (hereafter LOSVD) of the stars, it is necessary 
to correct the measured depths by a factor depending on 
the moments of the velocity distribution (Davies et al. 1993; 
Kuntschner 2000, 2004). The latter moments must be deter- 
mined using specialized code (Bender 1990; Saha & Williams 
1994; Pinkney et al. 2003; Merritt 1997; Kuijken & Merrifield 
1993; van der Marel & Franx 1993). These kinematical de- 
convolution routines have been used for some time and have 
undergone 2 major mutations. First, thanks to the increas- 



ing power of computers, it became affordable to swap back 
and forth from direct space to Fourier space, so that many 
disturbances such as border effects and saturation could be 
avoided. It became straightforward to mask problematic re- 
gions of the data, such as dead pixels, emission lines, etc... 
The second evolution of these codes allowed the use of mul- 
tiple superimposed stellar templates to best match the ob- 
served spectrum (Rix & White 1992; Cappellari & Emsellem 
2004). It has also been proposed to use single stellar popula- 
tions as synthetic templates, and this approach has proved 
to be useful in addressing the template mismatch problem 
(Falcon-Barroso et al. 2003). Moreover, this technique can 
save precious telescope time since it circumvents the need 
for observing template stars. 

In Ocvirk et al. 2005 (hereafter Paper I), we introduced 
STECMAP, a method for recovering non-parametrically the 
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stellar content of a given galaxy from its integrated light 
spectrum. Using STECMAP requires, as a preliminary, con- 
volving the data or models with the proper Point Spread 
Function (PSF), which can be of both physical {i.e. the 
stellar LOSVD) and instrumental (the instrument's PSF) 
origin. Adjusting the LOSVD to fit the data does not only 
constrain the kinematics of the observed galaxy but will also 
reduce the mismatch due to errors in the determination of 
the redshift or anomalous PSF, which is ultimately a neces- 
sary step when fitting galaxy spectra. 

Here we propose to constrain the velocity distribution 
simultaneously with the stellar content, by merging the kine- 
matic doconvolution and the stellar content reconstruction 
in one global Maximum A Posteriori likelihood inversion 
method. Hence, STECMAP becomes STECKMAP (STEUar 
Content and Kinematics via Maximum A Posteriori likeli- 
hood). In this respect, STECKMAP resembles the method 
proposed by e.g. Falcon-Barroso ct al. (2003), except that 
it takes advantage of the treatment of the stellar content by 
STECMAP. Together with the stellar age distribution and 
the age-metallicity relation, the LOSVD is described non- 
parametrically and the only a priori we use are smoothness 
and positivity. 

We also tentatively address the case of age-dependent 
kinematics, i.e. we try to recover the individual LOSVDs 
and ages of several superimposed kinematical subcompo- 
nents. This approach is motivated by the fact that galax- 
ies often display several kinematical conipoucnts. Ellipticals 
and dwarf ellipticals are for instance known to often harbor 
kinematically decoupled cores (De Rijcke et al. 2004; Balcells 
& Quinn 1990; Bender & Surma 1992), and spiral galaxies 
are usually assumed to be constituted of a thin and a thick 
disk, a bulge and a halo (Freeman & Bland-Hawthorn 2002). 
The variety of the dynamical properties of the components 
has a counterpart in their stellar content, as a signature of 
the formation and evolution of the galaxy. For instance, the 
halo of the Milky Way is believed to consist mainly of old, 
metal poor stars, while the bulge is more metal rich, and 
the thin disk is mainly younger than the bulge (Freeman & 
Bland-Hawthorn 2002). It is thus natural to let any stellar 
sub population have its own LOSVD. This possibility has 
been recently addressed by De Bruyne et al. (2004a,b), in 
a slightly different framework: they use individual stars as 
templates for the different components, while we propose 
to use synthetic SSP models. Such a method would allow 
us to separate the several kinematical components of galax- 
ies from integrated light spectra, and constrain for exam- 
ple their age- velocity dispersion and age-metallicity relation. 
The highly detailed stellar content and kinematical informa- 
tion that can be obtained for the Milky Way or for nearby 
galaxies that can be resolved into stars, such as M31 (Fer- 
guson et al. 2002; Ibata et al. 2004), could be extended to a 
larger sample of more distant galaxies. This technique could 
also be useful in detecting and characterizing major stellar 
streams in age and velocity from integral field spectroscopy 
of galaxies. 

In the whole paper we use the Pegase-HR SSP mod- 
els (Le Borgne et al. 2004) in order to illustrate and inves- 
tigate the behaviour of the problems through simulations 
and inversions of mock data. Indeed, Pegase-HR, with 
its high spectral resolution {R = 10 000), is an adequate 
choice for testing the recovery of detailed kinematical infor- 



mation in the form of non-parametric LOSVDs. The prob- 
lems and methods we describe arc however by no means 
specific to Pegase-HR (and its wavelength coverage) and 
STECKMAP could be used with any possible SSP model, 
depending on the type of data that is being analyzed. 

We will start with the modeling of the kinematics. Then, 
we will address the idealized linear problem of recovering the 
LOSVD when the stellar content is known, i.e. the template 
is assumed to be perfect. Section 4 deals with simultaneous 
age and LOSVD reconstruction of composite populations. 
Finally, section 5 investigates the case of age-dependent 
kinematics in a simplified context where the metallicity and 
extinction are known a priori. 



2 MODELS OF GALAXY SPECTRA 

In this section we present the modeling of galaxy spectra, 
taking into account the composite nature of the stellar pop- 
ulation, in age, metallicity and extinction, and finally its 
kinematics. 

2.1 The composite reddened population at rest 

We model the SED of the composite reddened population 
at rest using the ingredients defined in Paper I: 

J'rert(A) = Ut{E,\) jj^^^ mB{\t,Z{t))At, (1) 

where A(t) is the luminosity weighted stellar age distribu- 
tion, Z{t) is the age-metallicity relation, and B[\, t, Z) is the 
fiux-averaged single stellar population basis of an isochrone 
population of age t, /oxt the extinction law, and metallicity 
Z. We recall briefly the main properties of the Pegase-HR 
SSP basis we used in this paper. As mentioned earlier, spec- 
tral resolution is _R = 10 000 over the full optical domain 
AA = [4000, 6800] A, sampled in steps of 0.2 A. The models 
are available for metallicities Z £ [0.0001, 0.1] and consid- 
ered reliable between tmin = 10 Myr and tmax = 15 Gyr. 
The IMF used is described in Kroupa et al. (1993) and the 
stellar masses range from 0.1 M© to 120 Mq. The extinction 
law /ext was taken from Calzetti (2001). 

2.2 Model kinematics 

Stellar motions in galaxies define a LOSVD corresponding 
to projected local velocity distributions integrated along the 
line of sight and across one resolved spatial element. 

2.2.1 Global kinematics 

The motion of the stars can to first approximation be ac- 
counted for by assuming that the velocities of all stars of all 
ages along the line of sight axe taken from the same veloc- 
ity distribution (hence "global"). The model spectral energy 
distribution, (t>{X}, is the convolution of the assumed nor- 
malized LOSVD, (;(v), defined for v £ [vmin,Vmax] with the 
model spectrum at rest Frest(A). The convolved spectrum 
0(A) reads: 

^(A) = £7i^.e.(3^).(v)^, (2) 
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where c is the velocity of light. The above expression reads 
as a standard convolution 



/■Wmax ^ 

; / F{w — u)g{u)Au , 



with the following ropararnctcrization: 

w = ln(A) , it = hi(l + -), 

c 

F{W) = Fre.tCc"") =F,est(A), 

g{u) = g(c(c"-l))=fl(v), ^{w) = 

Wmin = ln(l + Xf^) , Uma:K = ln(l + 



(3) 

(4) 
(5) 

<^(A),(6) 
(7) 



2.2.2 Age-dependent kinematics 

Wo now allow the LOSVD to depend on the age of the stars. 
For simplicity, we consider here only unreddened mono- 
metallic populations, i.e. /ext(i^. A) = 1 and Z{t) = Zq. 
We introduce the age- velocity distribution, H(v,i), defined 
in [Vniin,Vmax] ^ [^min, tmax], which gives the contribution of 
stars of velocity and age in [v, v-|-dv] x [t,t + dt] to the total 
observed light. Thus, for a given age t, E(v, t) is the LOSVD 
of the single stellar population of age t. The age-velocity 
distribution H(v, t) is related to the stellar age distribution, 
A(^), by: 



H(v,f)dv = A(f), 



(8) 



The model spectrum of such a population thus reads: 

The above expression can be rewritten more conveniently 
^(w) = c/ / B{w -u,t)'S{u,t)dtdu, (10) 

•'Mmin •'tmin 



using the same reparameterization as in Sect. 2.2.1 and 

B{w,t) = B{e'",t,Zo) = B{\,t,Zo), (11) 

Eiu,t) s H(c(c"-l),t)=H(v,t). (12) 

In the rest of the paper, we will use exclusively the 
standard {i.e. reparameterized) convolutions as in Eq. (3) 
and Eq. (10). For readability, we will drop the superscript ~ 
and set the speed of light to one. 



3 KINEMATICAL DECONVOLUTION 

Sect. 2.2.1 shows that with proper reparameterization, the 
convolution of a model spectrum at rest, F{w), with the 
stellar LOSVD, g{u), reads as a standard convolution, given 
by Eq. (3). Finding the LOSVD when the observed spec- 
trum, (j>{w), and the template spectrum, F{w), are given is 
a classical deconvolution problem. Our goal here is not to 
discuss the respective qualities of the many different meth- 
ods available in the literature to solve this problem. Most 
rely on fitting the data while imposing some a priori on the 
LOSVD, I.e. they provide Maximum A Posteriori (MAP) es- 
timates of the LOSVD. Let us describe briefly our method 
to obtain such a solution with the purpose of coupling it in 
a later step with STECMAP. 



3.1 The convolution kernel 

Here we discretize Eq. (3) to obtain a matrix form defining 
the convolution kernel. We use an evenly spaced set 

{uj = Wmin + [j - 1/2) 6u;j = 1, 2, . . . ,p} 

spanning [umin,Wmax] with constant step 5u = ( 
'Umin)/p. We expand the LOSVD as a sum of p gate func- 
tions: 



9,0 



6u 



where 



e{x) 



1 if -1/2 < X ^ 1/2 
otherwise 



Injecting this expansion into Eq. (3) leads to: 



J29jFiw-Uj). 



(13) 



Similarly, we now sample along the wavelengths by integrat- 
ing over a small Sw. 



(14) 



where {wj-.j = 1,2, is a set of logarithmic wave- 

lengths spanning the spectral range with a constant step. 

Using matrix notation and accounting for data noise, 
the observed SED reads: 



Kg + e, 



(15) 



where y = {4'ij4'2, ■ ■ ■ ,<t'm) is the measured spectrum, 
and e = (ei, 62, . . . , Cm)^ accounts for modelling errors and 
noise. The vector of sought parameters g = {gi,g2, . . . , gp)~^ 
is the discretized LOSVD. The vector s = K • g is the model 
of the observed spectrum and the matrix K 

Kij = F{wi - Uj) , \/{i,j) e {1, . . . ,m} X {1, . . . ,p} , (16) 

is called the convolution kernel. 

The convolution theorem (Press 2002) states that the 
Fourier transform of the convolution of two functions is 
equal to the frequency-wise product of the individual Fourier 
transforms of the two functions. Applying this theorem 
yields another equivalent expression for the model spectrum 



s = JP""^ ■ diag(J!^ ■F)-T-s, 



(17) 



where is the discrete Fourier operator defined in Press 
(2002) as: 

= exp ( — {i - 1)U - 1)) , V(i,i) e [1, . . . mf ,(18) 



(19) 



Note that since m is the size of the template spectrum F, 
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the discretized LOSVD g, which is initially of size p needs to 
be symmetrically padded with zeros to the size m in order 
to transform the Toeplitz matrix into a circulant one. The 
diagonal matrix diag(jF ■ F) carries the coefficients of the 
Fourier transform of the model spectrum at rest F. This no- 
tation involving the Fourier operator, T, will be very useful 
for a number of algebraic derivations in the rest of the paper. 
In practice, from a computational point of view, it is more 
efficient to implement any forward or inverse Fourier trans- 
form through FFT. Similarly, the product diag(.7^ • F) • .7^ • g 
is in practice implemented as a frequency-wise product of 
the individual FFTs. 



3.2 Regularization and MAP 

A number of earlier publications have shown that the max- 
imum likelihood solution to Eq. (15) is very sensitive to the 
noise in the data e. Hence, in the spirit of Paper I, we choose 
to regularize the problem by requiring the LOSVD to be 
smooth. To do so, we use the quadratic penalization P(g) 
Eis defined by Eq. (29) in Paper I; 

P(g) =g^ L-g. (20) 

In the rest of the paper, the penalization is Laplacian, i.e. 
L = D2, where D2 is the discrete second order difference 
operator, as defined in Pichon et al. (2002). The objective 
function, Q^, to be minimized is given by: 

QM(g) =x'(y|g) + Mng)> (21) 

where the is defined by 

X'(y|g) = (y - s(g))^-W-(y - s(g)) . (22) 

The vector y is the observed spectrum and the weight ma- 
trix is the inverse of the covariance matrix of the noise: 
W = Cov(e)~^. The parameter 11 controls the smoothness 
of the LOSVD through its coefficients, g. It can be set on 
the basis of simulations (as described in Paper I) or auto- 
matically by GCV (Wahba 1990), according to the SNR of 
the data. In the latter case, the properties of the convolution 
kernel can be used to speed up the computation of the GCV 
function. Further regularization is provided by the require- 
ment of positivity, implemented through quadratic reparam- 
eterization. Minimizing yields the regularized solution 
g;j. Efficient minimization procedures require the analytical 
expression of the gradients of Q^, given in Sect. Al. 



3.3 Simulations 

We applied this deconvolution technique to mock data, cre- 
ated from Pegase-HR SSPs of several ages and metallici- 
ties, with i? = 10 000 at 4000-6800 A. In a first set of exper- 
iments, the model spectrum at rest was a solar metallicity 10 
Gyr single stellar population. It was convolved with various 
LOSVDs, both Gaussian and non Gaussian, with velocity 
dispersions ranging from 30 to 500 km/s. It was then per- 
turbed with Gaussian noise at levels ranging from SNR= 5 
to 100 per pixel, and deconvolved using the model spec- 
trum at rest as template {i.e. no template mismatch). In all 
cases, the LOSVDs are adequately recovered. Figure 1 shows 
the reconstruction of a Gaussian LOSVD, for SNR= 10 per 




Figure 1. Non parametric reconstruction of a Gaussian LOSVD 
for simulated data, CTi, = 100 km/s, SNR=10 per pixel. The model 
spectrum at rest is a 10 Gyr old solar metallicity single stellar pop- 
ulation with R=W 000 at 4000-6800 A. The template spectrum 
is identical, so that no template mismatch is allowed here. The 
thick line is the input model. The circles and the bars show respec- 
tively the median and the interquartiles of the recovered solutions 
for 500 realizations of the noise. 

pixel. However, there are necessarily some biases in the re- 
construction of the sharp features of the LOSVD. This is ex- 
pected since we introduced regularization via smoothing. To 
illustrate the relationship between regularization and bias, 
we performed a new set of similar simulations for a non- 
Gaussian LOSVD (sum of 2 Gaussians) with SNR=20 per 
pixel and varied the smoothing parameter /i. The results 
are shown in figure 2. The panels a and b correspond to 
^ = 10 while the panels c and d correspond to = 1000. 
The model, median and interquartiles of 500 reconstructions 
are displayed. We also plotted the whole set of 500 recov- 
ered solutions, in order to show the locus of the solutions. 
One can see that the biases of the median recontruction are 
reduced when lowering /x. The highest bump is correctly re- 
produced for /i = 10 while it is not for fj, = 1000. But on the 
other hand the solutions are much more widely spread when 
^ = 10. This means that most solutions taken from the set 
of low simulations can be very far from the model, while 
all the large fj, solutions lie reasonably close to the model. 

The regularization acts as a Wiener filter in the sense 
that it damps the high frequency components of the solution. 
Regularization improves the significance of an individual re- 
construction (it will nearly always lie reasonably close to the 
model), at the cost of introducing a bias. 

3.4 Age and metallicity mismatch 

We take advantage of the large range of ages and metallic- 
ities of single stellar populations covered by Pegase-HR 
to shortly illustrate the effects of template mismatch on 
LOSVD determinations. In this section we show the results 
of a large number of simulations aiming at characterizing the 
error made when a wrong template is chosen for the kine- 
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(a)SNR=20 mu=10 (b) SNR=20 mu=10 
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Figure 2. Impact of the smoothing parameter. Reconstruction of a non-Gaussian LOSVD for simulated data with SNR=20 per pixel. 
Top: = 10, bottom: fi = 1000. Left: the thick line is the input model LOSVD. The circles and the bars show respectively the median and 
interquartiles of the reconstructed LOSVDs for 500 realizations of the experiment. Right: the whole set of 500 solutions is displayed, with 
the model as a thick white line, in order to give the reader a sense of what individual solutions look like. The figures show the tradeoff 
between bias and reliability of the reconstruction: for small /i the median reconstruction is unbiased but the individual reconstructions 
are very noisy. For large /i the median reconstruction is slightly biased but all the reconstructions are reasonable. Hence, the significance 
of an individual reconstruction is improved by regularization at the cost of introducing a bias. 



matical inversion of data. For this purpose, mock data were 
created by convolving a single stellar population of age, ao, 
and metallicity, Zq, with a centered Gaussian LOSVD of dis- 
persion (T„ = 50 km/s. It was perturbed by Gaussian noise 
corresponding to SNR= 100 per pixel and then deconvolved, 
using as template a single stellar population of age ai, and 
metallicity Zi . The spectral resolution and wavelength range 
are the same as in Sect. 3.3. Fig. 3 shows the error on the 
measured velocity dispersion. The latter is measured as the 
r.m.s of the reconstructed LOSVD. If the parameters of the 
template are different from those of the model, the velocity 
dispersion error increases very quickly. The age metallicity 
degeneracy is visible as a valley of smaller error, following 
the upper-left to bottom-right diagonal of the figures. Of 
course, the distance between the model and the mock 



data follows a similar 2D distribution, and will lead to the 
rejection of highly mismatched LOSVD estimates. However, 
in practice, it is usually not straightforward to quantify all 
the sources of error. It is thus somewhat arbitrary to set an 
upper limit of for the admissible solutions, and the error 
on the kinematics is thus hard to quantify. This experiment 
illustrates in this context the long known issue that when 
the correct model is not available, large errors on the de- 
termination of kinematics are expected. In order to reduce 
the error in the estimates of the kinematical properties of a 
stellar assembly, it is necessary to allow for a wide range of 
modulations of the template. This is naturally achieved by 
making the non parametric stellar content account for the 
changes of the template, as discussed in the next section. 
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Figure 3. Velocity dispersion error as a function of the age and rnctallicity of the template single stellar population. Contours show 
regions of increasing velocity dispersion error. In each experiment, the age and metallicity of the original model template is shown as a 
thick cross, and the model LOSVD is a Gaussian with zero mean and 50 km/s dispersion. The velocity dispersion error is minimum when 
the template's age and metallicity are similar to the model's. The error increases quickly when the template's parameters differ from 
the model's, also in the age-metallicity degenerax;y direction (upper left-bottom right diagonal). It increases even faster in the direction 
orthogonal to the age-metallicity degeneracy. 



4 RECOVERING STELLAR CONTENT AND 
GLOBAL KINEIVIATICS 

The mixed inversion described in tliis section couples tlie 
recovery of both the stellar content and the kinematics, 
thereby turning STECMAP into STECKMAP. Proper ap- 
plication of this method provides an interpretation of the 
observed object in terms of stellar content and kinematics. 

4.1 Inverse problem 

For a given model spectrum at rest, -Frcst(A), and a LOSVD, 
g{v), the emitted SED, (j>{X) is given by Eq. (2). We now 
wish to account also for the additional variables involved in 
-Frcst, given by Eq. (1), namely the stellar age distribution, 
A(t), the age-metallicity relation Z{t), and the color excess 
E{B — V) = E. Injecting Eq. (1) into the convolution Eq. (3) 
yields the emitted SED: 

ct>{w) = jj /ext [E, w - u)A{t)B{w - u, t, Z{t))g{u) At Au , 

(23) 

Solving Eq. (23) for A, E and g when i^, /ext and B are 
given is the inverse problem we are tackling here. 

4.2 Discretization and parameters 

Expanding the two time-dependent unknowns A(t) and Z(t) 
as a sum of n gate functions and injecting into Eq. (1) yields 
the discrete model spectrum at rest: 

F = diag(fext(S))-B-x, (24) 

This discretization is explained in details in Sec. 5 of Paper 
I. Similarly, we develop the LOSVD g{u) as a sum of p gate 
funtions as in Sect. 3. Note that the reddened model at rest 
plays the role of the stellar template in a classical kinematic 
convolution. Injecting Eq. (24) into Eq. (17) thus allows us 
to express the model spectrum, s, as 

5 = 0^ ■ diag(.F • diag(fext (-E)) • B • x) • • g , (25) 



However, here, the template is this time modulated by the 
unknowns describing the stellar content. 

4.3 Smoothness and metallicity constraints 

The discrete problem of finding the stellar age distribution 
X, the age-metallicity relation Z, the extinction E and the 
LOSVD g for an observed spectrum y and given an extinc- 
tion law /oxt and a SSP basis B is of course likely to be very 
ill-conditioned since it arises as the combination of several 
ill-conditioned problems. It therefore requires regularization. 
We also want the metallicity of the components to remain 
in the model range. We use the standard penalization P 
and the binding function C defined in Paper I to build the 
penalization for this problem: 

P^ = /ixP(x) + nzP{Z) + McC(Z) + n^P{g) , (26) 

where fi = (^x, Mz, PC A*")- Again, we choose L = D2 as 
defined in Pichon et al. (2002), so that the penalization P 
is actually Laplacian. The objective function, Q^, is now 
defined as: 

=x'(s(x,Z,S,g))+P^(x,Z,i5,g). (27) 

and its partial derivatives are given in Sect. A2. Note that 
there is in principle an additionnal formal degeneracy for 
this inverse problem. If the set (x, Z,i?,g) is a solution to 
(23), then (ax, Z, _E, g/a) is also a solution for any scalar 
Q, because the age distribution x and the LOSVD g are 
not explicitly normalized in this formulation. However, the 
adopted regularization lifts this degeneracy. The penaliza- 
tion function P is quadratic (P(qx) = a'^P(x)). Thus, if x 
or g is too large in norm, the solution is unattractive. Prac- 
tically, the algorithm reaches a solution where x and g are 
similar in norm. In any case, this degeneracy would easily be 
remedied by adding a normalizing term to the penalization 
Pp of the form ||x|| — 1, which would force the discretized 
stellar age distribution x to have unitary norm. Following 
the same principle, one could equivalently choose to normal- 
ize the LOSVD rather than the stellar age distribution. 
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4.4 Simulations 

Let us now test the behaviour of STECKMAP by applying it 
to mock data. The latter were produced using an arbitrary 
stellar age distribution x, an age-metallicity relation Z, a 
LOSVD g and an extinction parameter E. Several simula- 
tions were performed with various input models: bumpy age 
distributions, increasing or decreasing age-metallicity rela- 
tion and extinctions, Gaussian and non Gaussian wide or 
narrow LOSVDs, in various pseudo-observational contexts. 
Figure 4 shows the results of two of these experiments. In 
the top line, the model is a young metal-poor population su- 
perimposed to an older metal-rich population. In the bottom 
panels, the model has a rather constant stellar age distribu- 
tion, a non-monotonic age-metallicity relation and a strongly 
non Gaussian LOSVD. In both cases the 3 unknowns are cor- 
rectly recovered. In these examples, the data quality mimics 
that of the best Sloan Digital Sky Survey galaxies: the reso- 
lution is Ji « 2000 and SNR= 30 per ?»1 A pixel. The wave- 
length domain of Pegase-HR is however narrower than the 
SDSS's. These simulations simply aim at demonstrating the 
generally good behaviour of the method, and show that ac- 
counting for the kinematics does not fundamentally weaken 
the constraints on the stellar content. For a more thorough 
study of the informational content of the Pegase-HR wave- 
length range, the reader can refer to the systematic double 
burst simulations with variable spectral resolution and SNR 
per A performed in Paper I. 



5 RECOVERY OF AGE-DEPENDENT 
KINEMATICS 

In this section we present an implementation of the recov- 
ery of age-dependent kinematics, i.c the situation when each 
sub-population has its own LOSVD. In this experiment, we 
restrict ourselves to the case where the stellar populations 
have a known metallicity and are seen without extinction. 
This choice is mainly motivated by the numerical cost of 
such a large inversion procedure. The modeling is given by 
Eq. (10). Finding the age- velocity distribution H(u, i) when 
the mono-metallic basis B and the observed spectrum 4> are 
given is the inverse problem. It arises as the combination of 
a linear age inversion and a kinematical deconvolution. 



5.1 A sum of convolutions 

The age- velocity distribution, H(u, t), is expanded as a linear 
combination of normalized 2D gate functions 9ij{u,t): 



St 



In other words, H(u, t) is represented by a 2D array v of size 
(p, n) , i.e. p is the size of each LOSVD and n is the rmnibcr 
of age bins. The linear step in u is 5u and the step in t is St. 
By injecting the expansion into Eq. (10) we obtain 

4i{'w) = I I VijSij {u, t)B{w — u, t)dtdu , 



i=i j=i 

p n 



i=i j=i 



(28) 



As in the previous sections, Bj (u) is a time-averaged single 
stellar population of age ± | St. We then discretize along 
wavelengths by averaging over small Sw. 



(t>k 



^ P n J. 



i:=l j = l 
p n 



^^VijBj{wk-Ui), 

i=l 3 = 1 



(29) 



where {wj)je{o,...,m} is a set of constant step logarithmic 
wavelengths. The above expression also reads in matrix form 
as a sum of kernel convolutions. Finally, the model SED of 
the emitted light reads: 



where s = {4>i, 



^), Vj = {Vlj 
K22j 



(30) 



mlj 



m2j 



V2j, 

K 



, Vpj ) and 



2pj 



with 



Kikj = Bj{wk — Ui 



(31) 



(32) 



With this notation, Kj and Vj are respectively the convolu- 
tion kernel and the LOSVD of the sub-population of ago tj , 
and the model spectrum y is the sum of the convolution of 
the kernel of each sub-population by its own LOSVD. 

5.2 2D age-velocity smoothness constraints 

In the previous sections, the unknowns were mono- 
dimensional functions of time or velocity. Here, the unknown 
is a 2D distribution, and we thus have to implement a 2D 
smoothing constraint. We wish to allow the smoothness in 
age to be distinct from the smoothness in velocity. We thus 
construct two penalizing functions. Pa and Pv , relying on the 
standard function P. Pa computes the sum of the Laplacians 
of the columns of v while P„ computes the sum of the Lapla- 
cians of the lines of v. The smoothness in the direction of 
the velocities (respectively ages) is set by /x„ (respectively 
/ia). We define the vectors Vj — {vij,V2j,--- ,Vpj) as the 
columns of v, i.e. the LOSVD s of the subpopulations. We 
similarly define the v' = (wa, «i2, • ■ ■ , Vin) as the lines of v. 
With this notation, the penalization reads: 

P^(v) = HaPa{v) + l.lvPv{-v) , 
p n 
= HaY,P{v')+l.ivY.P{Yj). (33) 

i=l j=l 

The objective function, Q^, is now fully specified as = 
+ P/j.- Its gradients are given in appendix A3. We choose 
here the smoothing parameters, n = {na,l^v), on the basis 
of simulations. 



5.3 Simulations of a bulge-disk system 

We studied the feasibility of separating two age-dynamically 
distinct populations, i.e. two components which do not over- 
lap in an age-velocity distribution diagram, in a regime of 
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Figure 4. Reconstruction of the stellar age distribution, age-metallicity relation and LOSVD for simulated SDSS-likc data with SNR=30 
per pixel. The thick histogram is the input model. The circles and the bars show respectively the median and the interquartiles of the 
recovered solutions for 50 realizations. 



very high quality model and data. We performed simulations 
in the idealized case of a very simplified spiral galaxy consist- 
ing of a bulge-disk system of solar metallicity seen without 
extinction at some intermediate inclination, in two observa- 
tional contexts. The corresponding ages and projected kine- 
matical parameters are given in Table 1. The resolution of 
the pseudo-data is = 10 000 at 4000 - 6800 A, and the 
SNR is 100 per 0.2 A pixel. 



• Case 1: The galaxy is resolved, and the fiber aperture 
is small compared to the angular size of the galaxy. The line 
of sight is offset by a couple kpc from the center along the 
major axis. The projected model age-velocity distribution 
involves 2 superimposed components: an old, non-rotating 
kinematically hot population representing the bulge, and a 
young, rotating, kinematically cold component. The model 
and the median of 30 reconstructions are shown in Fig. 5. 
The separation of the components is clear and their param- 
eters can be recovered with good accuracy, considering the 
difficulty of the task. 

• Case 2: The galaxy is unresolved. The difi'erence with 
the former situation is that because of the spatial inte- 
gration, both age-velocity distributions are centered. For a 
given dynamical model, the projected dispersion of the disk 
component depends on its inclination. Fig. 6 shows that the 
separation is successful and that the ages and integrated 
kinematical properties of both components can be measured. 



Case 1 



ye(km/s) I i7v(km/s) | age (Gyr) 



Bulge 





100 


8 


Disk 


120 


30 


0.5 


Case 2 




Bulge 





150 


8 


Disk 





50 


0.5 



Table 1. Projected kinematical parameters and age of the model 
bulge-disk system used to produce the simulations of Fig. 5 and 
Fig. 6. Vc (respectively CTv) is the rotation velocity (respectively 
the velocity dispersion) projected on the line of sight. 



6 CONCLUSIONS 

The non-parametric kinematical deconvolution of a galaxy 
spectrum is efficiently performed using a MAP formalism 
(Sect. 3). Regularization through smoothness requirements 
and positivity improve significantly the behaviour of the in- 
version with respect to noise in the data. This improvement 
occurs at the cost of introducing some bias in the recon- 
structed LOSVD, but this bias remains reasonable. Strong 
non-Gaussianities of LOSVDs are reliably detected from 
mock data generated using Pegase-HR SSPs for SNR down 
to 20 per 0.2 A pixel. 

When the template does not exactly match the model 
spectrum at rest, i.e. there is some template mismatch, 
the error on the velocity dispersion increases very quickly 
(Sect. 3.4). For example, in our experiments, where av = 50 
km/s with 7? = 10 000 data, the error on the measured veloc- 
ity dispersion amounts up to 10-20 per cent if the template 
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model median reconstruction 




v[km/s] v[km/s] 

Figure 5. Model (left) and median reconstruction (right, fx 30 realizations) of a stellar age-velocity distribution from SNR=100 per 
0.2 A pixel pseudo-data at 4000—6800 A. The model stellar age-velocity distribution mimics that of a simplified spiral galaxy seen with 
intermediate inclination. The old broad component can account for the bulge population while the young narrow rotating component 
represents a thin disk population. The projected kinematical parameters of the model are given in Table 1 (case 1). The different 
kinematical components are well separated and clearly identifiable. 

model median reconstruction 




v[km/s] v[km/s] 

Figure 6. Same as Fig. 5 but for an unresolved simplified spiral galaxy with projected and spatially integrated kinematical parameters 
given in Table 1 (case 2). The velocity dispersion of the integrated young disc component depends on the inclination angle. The bulge 
and the disc are are well separated and clearly identifiable. Their respective velocity dispersions and ages are reliably recovered. 



diflters from the model by more than 0.3 dex in age and 
metallicity, perpendicular to the age-metallicity degeneracy. 

The formal similarity between the non-parametric kine- 
matical deconvolution and the recovery of the stellar content 
allows us to merge both processes in a "mixed" inversion 
where the observed spectrum is fitted by determining the 
stellar content and the kinematics simultaneously (Sect. 4). 
This circumvents the need for iterations where kinematical 
and stellar content analyses are carried out one after the 
other, until convergence is reached; this provides an efficient 
method to analyze large sets of data. 

Satisfactory reconstructions of the stellar age distribu- 



tion, the age-metallicity relation, the extinction and the 
global LOSVD were obtained from mock data down to 
R = 2000, SNR= 30 per 1 A pixel in the 4 000-6 800 A range 
(simulating SDSS data in the Pegase-HR range) , indicating 
the good behaviour of the method. Since in our simulations, 
the introduction of the kinematics into STECMAP did not 
affect the recovery of the stellar content, we consider that 
the error estimates and separability analysis given in Paper 
I remain valid. 

In a more exploratory part of this work, we showed the 
feasibility of recovering age-dependent kinematics in a sim- 
plified mono-metallic unreddened context (Sect. 5). We were 
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able to separate the bulge and disc component of a simpli- 
fied model spiral galaxy in integrated light provided very 
high quality data (SNR=100 per 0.2 A pixel in the optical 
domain) and models are available, i.e. we constrain both 
components in velocity dispersion and ago. This separation 
was also carried out successfully in the setup corresponding 
to an unresolved galaxy. 

Farther investigations are needed to extend this tech- 
nique to a regime where the metallicity and extinction are 
unknown. Wc expect that letting the metallicity be a free 
parameter would certainly lead to a more degenerate prob- 
lem, as shows the degradation of the resolution in age found 
in Paper I compared to fixed metallicity problems. On the 
contrary, wc do not expect the addition of the extinction 
as a free parameter or a more complex form of extiction 
law or flux calibration correction, possibly non parametric, 
to deteriorate the conditioning of the problem. The results 
arc encouraging, and the feasibility of such age-dependent 
kinematics reconstructions deserves to be tackled in realistic 
specific pseudo-observational regimes in the future. 

As mentioned in Paper I, the SSP models were consid- 
ered to be perfect and noiseless. It still has to be investigated 
how instrumental error sources such as flux and wavelength 
calibration error, additive noise, contamination by adjacent 
objects, and, equally important, model errors, can affect the 
robustness of such sophisticated interpretations. 



PERSPECTIVES 

STECKMAP wUl be very useful to interpret data of large 
spectroscopic surveys, complete or in progress, such as 2DF- 
GRS\ SDSS^ DEEP2-\ or VVDS'', especially where both 
constraints on the stellar content and the dynamics are re- 
quired. STECKMAP's analysis of the spectroscopic survey 
data or of a SNR selected subsample, combined with sur- 
vey photometry could provide estimates of the stellar and 
dynamical masses (which must be corrected for flber aper- 
ture though), thereby allowing astronomers the prospect of 
investigating the dark matter content in galaxies on a star 
tistically significant sample, in the spirit of Padmanabhan 
et al. (2004). 

The application of agc-dopondent kinematics to integral 
field spectroscopy data from, for example SAURON (Bacon 
et al. 2001; de Zeeuw et al. 2002), OASIS (McDermid et al. 
2004), MUSE (Henault et al. 2003) or MPFS (Chilingarian 
et al. 2004) could significantly boost the amount of informa- 
tion extracted from this data. 

The inner paxts of elliptical or dwarf elliptical galaxies 
have shown via adaptive optics new kinematically decoupled 
structures (cores or central disks), which were precedently 
unresolved (McDermid et al. 2004; Bacon et al. 2001). Sim- 
ilarly, if decoupled structures are unresolved and remain so, 
even with adaptive optics, it may still be possible to separate 
components in age-velocity space. Hence, the technique pre- 
sented in Sect. 5 extends the range of investigation for the 
inner components of galcixies even further in redshift and 

^ http://www.mso.anu.edu.au/2dFGRS/ 

^ http://www.sdss.org/ 

^ http://www.deep.berkeley.edu/ 

^ http://www.oamp.fr/virmos/wds.htm 



distance with the current generation of instruments. The 
faint, generalized counterparts of kinematically decoupled 
cores, i.e. stellar streams generated by minor merging and 
accretion of satellites, may also be detectable by an age- 
dependent kinematics reconstruction in systems which can 
not be resolved into stars, provided that they are sufficiently 
distinct from the bulk stars of the galaxy in the age- velocity 
space. This will enlarge the sample of galaxies for which such 
detailed information is available, and may make it statisti- 
cally significant. 
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APPENDIX A: GRADIENT COMPUTATIONS 
Al Kinematic deconvolution 

In this section we derive the gradient of with respect to 
the LOSVD g. First, we rewrite the term as: 

X^ = r^-W-v, (Al) 

where the residuals vector r is defined by 

r = y-JP-'-diag(^-F).:r.g, (A2) 

The derivative of the then reads: 

= -2JF* . diag(:F . F)* • :F . W • r , (A3) 

where the asterisk * denotes the complex conjugate. Since 
the stellar template and the LOSVD can play symmetrical 
roles in Eq. (17) we can also write the derivative of x^ rela- 
tively to the stellar template: 

^ = -2.^*-diag(.^-g)*-.F-W.r, (A4) 

This expression will be useful for later derivations of gradi- 
ents for more complex problems in the following appendices. 



A2 Gradients of the mixed inversion 

Here we show how to obtain the partial derivatives of 
Qfj, = + Pfj. defined in Sect. 4. Given that writing 
the derivatives of the penalizing functions is straightfor- 
ward, we will in this appendix focus on the gradients of the 
X^. In the mixed inversion, the reddened model spectrum 
at rest plays the role of the stellar template F in the clas- 
sical kinematic deconvolution of equation (15). dx^ /dg can 
thus be obtained by replacing F <— diag(fext(£^)) • B • x in 
equation (A3). 



dg 



= -2JF*.diag(:F.diag(fext(-B))-B • x)*-:F.W • r , (A5) 



where r = y — s is the residuals vector, with s as given by 
Eq. (25). To obtain the other partial derivatives we use the 
following relation. For any parameter a we have 



da 



OF 



da 



(A6) 



The first term dx^/dF is given by Eq. (A4) while the second 
term reads, considering each unknown: 

dF 



diag(fext) - A, 



II = diag(x) • ^ • diag(fext) , 

OF ,5fext^ . 

dE = d>ag(^)-A.x, 

with the same notation as in the appendix of the STECMAP 
paper. 



(A7) 
(A8) 
(A9) 



A3 Gradients for the age-dependent kinematics 
recovery 

Again, we focus on the partial derivatives of the x^- Us- 
ing Eq. (17), the model can be rewritten using the Fourier 
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operator 



n 



S = 



^^*-diag(^-B,).^.v, 



J 1 



(AlO) 



3 = 1 



where Bj is the discretized time-averaged single stellar pop- 
ulation of age ltj-i,tj], The derivatives of relatively to 
V can be derived directly from Eq. (A3) since the model is 
just a sum of convolutions. Replacing F •<— Bj and g •<— Vj 
yields the gradient of the x^'- 



^ = -2JF*.diag(j^-F)"'-:r.WT, (All) 



with the residuals vector r = y — s. Finally, the derivative 
of relatively to v is the matrix defined by 



dv 



( 



(A12) 



